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O ■ Abstract 

(N' 

►^ ■ We propose a general method for measuring an arbitrary observable of a 

lO ■ multimode electromagnetic field using homodyne detection with a single local 

z^ \ oscillator. In this method the local oscillator scans over all possible linear 

C^ \ combinations of the modes. The case of two modes is analyzed in detail and 

the feasibility of the measurement is studied on the basis of Monte-Carlo 
Q>^ ■ simulations. We also provide an application of this method in tomographic 

testing of the GHZ state. 
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I. INTRODUCTION 

Optical homodyne tomography is a well-established quantitative method for measuring 
the quantum state of radiation and for obtaining the expectation value of arbitrary ob- 
servables of the field ||l]-§ (for a review see Ref. |Q). The success of optical homodjTie 
tomography has stimulated research relating to the state-reconstruction procedures in other 
fields, such as in the realm of atomic 0, molecular 0, and ion-trap 0] physics. As a matter 
of fact, the tomographic method is a kind of universal detection technique [Q, with which 
one can measure any observable O of the field by averaging a suitable unbiased estimator 
£^[O](x,0) over the homodyne data x at random phase values 0. Single-mode homodyne 
tomography can be immediately generalized to multimode fields. For factorized multimode 
operators O = Oi ® O2 ® ■ ■ ■ ® On the corresponding estimator is just the product of the 
estimators for each of the single-mode operators Oi, Oi, . . . , On- By linearity the estimator 
can then be extended to generic multimode operators. However, such a simple generaliza- 
tion requires a separate homodyne measurement for each of the modes, which cannot be 
achieved in practice when the modes of the field are not spatio-temporally separated. For 
this reason, tomographic methods have been devised which either use only a single local 
oscillator (LO) [^, or avoid the use of conventional homodyne detection |T^. However, both 



the methods work for only two modes of the field, and the self- homodyne method of Ref. |T^ 
is suitable only in special experimental situations (e.g., in the tomography of parametrically 
downconverted radiation). Therefore, a more general multimode tomographic method is 
needed, especially in consideration of the possibility of a precise analysis for pulsed fields, 
for which the problem of mode matching between the LO and the detected fields (deter- 
mined by their relative spatio-temporal overlap) |jll[] gives a detrimental contribution to the 
overall quantum efficiency. 

In this paper we propose a general method for measuring an arbitrary observable of 
the multimode electromagnetic field, which uses homodyne detection with a single LO. We 
provide the rule for evaluating the "unbiased estimator" of a generic multimode operator. 
The quantum expectation value of the operator can then be obtained for any unknown state 
of the radiation field through an average of this estimator over the homodyne outcomes 
that are collected using a single LO which scans over different linear combinations of the 
incident modes. The paper is organized as follows: In Sec. II we present the general method 
for obtaining the estimator pertaining to an arbitrary multimode operator. Upon averaging 
this estimator over the homodyne outcomes, one obtains the quantum expectation value 
of the corresponding operator. We specialize to observables corresponding to the matrix 
elements of the multimode density operator and to the total number of photons. In the two- 
mode case we explicitly derive the estimator for the four-dimensional Q-function and for 
the moments generating function of the total number of photons. In Sec. HI we investigate 
the experimental conditions for extracting the joint photon-number probability and the 
distribution of the total number of photons for two-mode quantum states. We present 
the results of some Monte-Carlo simulations for the twin-beam state that is produced by 
nondegenerate parametric amplification (spontaneous downconversion) . We average the 
estimators obtained in Sec. II over the homodyne data that are distributed according to 
the theoretical homodyne probability evaluated in the Appendix. The simulations show 
that the measurement is feasible for quantum efficiency values of the homodyne detector 



in the 80-90% range and with the number of experimental data samples of order 10^-10^. 
In Sec. IV we show an application of our method in measurement of the three-particle 
maximally-entangled state called the GHZ state. In such a case the number of radiation 
modes is six and a more suitable arrangement of the tomographic machine requires the use 
of three LO's. The results of Monte-Carlo simulations show that for homodyne detectors 
with quantum efficiency value r] = 85% one needs about 10^ data samples to reconstruct the 
state with a relatively small statistical error. Finally, some conclusions are drawn in Sec. V. 

II. THE GENERAL METHOD 

For a single-mode radiation field one has the following resolution of the identity on the 
Hilbert- Schmidt space: 
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where O is a Hilbert-Schmidt operator, the integral is extended to the complex plane C for 
z, and D{z) = exp{za^ — z*a) denotes the displacement operator for the field mode with 
annihilation and creation operators a and a"!", respectively, having the commutation relation 
[a, a^ = 1. Equation (|l]) simply follows from the orthogonality relation for displacement 
operators Tt[D{z)D''{z')] = 62{z—z'), S2{z) denoting the Dirac delta-function on the complex 
plane. Equation (|l]) is the starting point of our method; it can be easily generalized to any 
number of modes as follows: 
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where a; and aj, with / = 0, . . . , M and [«;,«[/] = Sw, are the annihilation and creation 
operators, respectively, of the M + 1 independent modes, and O now denotes an operator 
over all the modes. Using the following hyper-spherical parameterization for zi & C 

zq = -kuo{e)e''^° = -fce^'^^cos^i , 

z^ = -kui{d)e''f'' = -ke^'f'^ sin 01 cos 02 , 
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ZM-i = ^k UM-i{9)e'^''-' = ^k e'^''-' sin 9i sin 62... sin 9m-i cos 9m 

zm = ^k UM{9)e''^^' = ^k e'"^'' sin 9i sin ^2 • • • sin 9m-i sin 9m , 

where k G [0,oo); ipi G [0,27r] for / = 0, 1, . . . , M; and 9i G [0, 7r/2] for / = 
Eq. (^ can be rewritten as follows: 
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O = j dfi[4j] J dfi[e] j^ dk M — Tr{OexpHfcX(^,^)]}exp[zA;X(^,V^)] . (4) 



Here we have used the notation 
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Notice that, thanks to the parameterization in Eq. (^, where X]i^'"f(^) — ^i ^^^ ^^^ 
the commutation relation [A{6 , ifj) , A^ {6 , tp)] = 1, which imphes that A{9,ip) and A^{9,ip) 
themselves are annihilation and creation operators, respectively, of a bosonic mode. Also, 
by scanning all values of 6i G [0, tt/2] and ipi G [0, 2tt], all possible linear combinations of the 
modes described by annihilation operators a;, with / = 0, . . . , M, are obtained. 

For a single mode of the radiation field the experimental homodyne probability distri- 
bution of a field quadrature with quantum efficiency r^ < 1 is a Gaussian convolution with 
variance A^ = ^j^ of the ideal probability distribution. Therefore, for the quadrature oper- 
ator X{6,ip) in Eq. (^, one has the following identity for the moments generating function 

(e^^^) = exp (^—^kA f dxe'''''pr,{x;e,^) , (8) 

— * — * — * — * 

where Pn{x]9,ip) denotes the homodyne probability distribution of the quadrature X{6,il)) 
with quantum efficiency r]. Generally, r/ can depend on the mode itself, i.e., it is a function 
7] = rjlOyip) of the selected mode. In the following, for simplicity, we assume r] to be mode 
independent, however. By taking the ensemble average on each side of Eq. (^ and using 
Eq. (W) one has 
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where, for a given operator O, the function Sn[0] (x; 6*, ip) of x, 6*, ip has the following analytic 
expression 

^ ^ ^M+l /•+00 ^ ^ 
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with 



K=^^. (11) 



Equation ([T0| ) is the central result of this paper. For any given operator O it provides the 
"unbiased estimator" to be averaged over all homodyne outcomes of the quadrature X{6,ijj) 
in order to obtain the ensemble average (O) for any unknown state of the radiation field. 
The homodyne outcomes for X{6, ip) can be obtained by using a single LO that is prepared 
in the multimode coherent state ^fi^lji) with 7/ = e^'^'ui{6)K/2 and K ^ 1. In fact, in 
this case the rescaled zero-frequency photocurrent at the output of a balanced homodjTie 
detector is given by 
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which corresponds to the operator X{6,ip). In the limit of a strong LO {K -^ 00), all 
moments of the current i correspond to the moments of X(^, ip), and the exact measurement 
of X{6, tp) is then realized. Notice that for modes a; with different frequencies, in the d.c. 
photocurrent in Eq. ( [T^ ) each LO with amplitude 7^ selects the mode ai at the same 
frequency (and polarization). For the effect of less-than- unity quantum efficiency, previous 
considerations on Eq. (|^) apply. 

In order to obtain the ensemble average in Eq. (^ as an experimental average over 
the homodyne outcomes, one needs to satisfy the validity conditions of the central-limit 
theorem. Since in the strong LO approximation the probability pjj{x; 6, ip) must decay as a 
Gaussian for large x, it follows that the integral in Eq. (P) can be experimentally sampled 
for any a priori unknown probability distribution Pri{x; 9, if)) if Sn[0] {x; 9, ip) increases slower 
than exp(/€a;^) for large x, and is bounded for |x| < +00. In this case one is guaranteed 
that the integral in Eq. (^) can be statistically sampled over a sufficiently large set of 
data. The average values for different experiments will be Gaussian distributed around the 
mathematical expectation in Eq. (^, allowing estimation of the confidence intervals, which 
will decrease as the inverse square-root of the number of experimental data. In general, 
the boundedness of £r\0]{x]9,ip) for |x| < +00 will pose lower bounds for the quantum 
efficiency rj below which the measurement cannot be performed, similarly to what happens 
in the one-mode case []T2| . This limitation is due to the generality of the method, which is 



perfectly unbiased, and makes no a priori assumption on the state of the radiation field, 
the only approximation being that of a strong LO. This should be contrasted with other 



methods, as the maximum entropy method |]T3[ or the maximum likelihood method |I^,p!5| 



which do not suffer such limitation on the quantum efficiency; however they are generally 
biased and based on assumptions for the state of the radiation field. 

Equation (P) can be specialized to some observables O of interest. In particular, one can 
estimate the matrix element {{ni}\R\{mi]) of the multimode density operator R. This will 
be obtained by averaging the following estimator: 
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where fii = max(m;,n/), z// = m.m{mi,ni), and L'^{z) denotes the customary generahzed 
Laguerre polynomial of variable z. For diagonal matrix elements, Eq. (13) simplifies to 
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with Ln{z) denoting the customary Laguerre polynomial in z. Notice that the estimator in 
Eq. (|I^) does not depend on the phases iph only the knowledge of the angles 9i is needed. 
Using the following identity for the Laguerre polynomials \W\: 
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from Eq. (|Tj) one can easily derive the estimator for the probability distribution of the total 
number of photons A^ = X]i=o '^I'^i 
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where \n) denotes the eigenvector of A^ with eigenvalue n. Notice that the estimator in 
Eq. (^) does not depend on any of the phases ifji or the angles 6*^, and thus their knowledge 
is not needed in the measurement of the probability distribution of the total number of 
photons. 

Now we specialize to the case of only two modes (i.e., M=l and 6' is a scalar 9). The joint 
photon-number probability distribution is obtained by averaging the following estimator: 
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Analogously, using Eq. (|TOD one can derive the following estimator for the four dimensional 
Q-function: 
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where \a, P) with a, /3 G C denotes a two- mode coherent state, and $(a, b; z) is the customary 
confiuent hypergeometric function of z. The estimator (|T6|) for the probability distribution 
of the total number of photons can be written as 
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For the total number of photons one can also derive the estimator for the moment generating 
function, using the generating function for the Laguerre polynomials |T6|. One obtains 
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where we have denoted by a and h the annihilation operators of the two modes. For the first 
two moments one obtains the simple expressions 

£r,[o)a + 6^6] (x; 6, tfjo, tpi) = Ax^ + - ~ 2 , (21) 

£J{a^a + b^bf]{x; 9, V^o, V^i) = 8x^ + f — - 20 ) x" + -^ - — + 4 . (22) 
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It is worth noting that analogous estimators for the photon-number difference between the 



two modes are singular and one needs a cutoff procedure, similar to the one used in Ref. |T^ 
for recovering the correlation between the modes by means of the customary two-mode 
tomography. The singular behavior of the estimators for the photon-difference operators 
can be understood simply from the fact that to extract information pertaining to a single 
mode only one needs a delta-function at ^ = for mode a, or 6 = 7r/2 for mode b, and, in 
this case, one could better use the original one-mode tomography method [|12| by setting 
the LO to the proper mode of interest. 

Finally, we note that for the case of two-mode tomography the estimator S^j can be 
averaged by taking the integral 

^"d^o /"^"dV^i [^ d(cos20) '•+°° 



(0)= ^ ^ n dxp,ix■,e,^Po,^Pl)£,[0]ix;e,^o,A) (23) 
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over the random parameters cos(26'), tpQ, and tpi. For example, in the case of two radiation 
modes having the same frequency but orthogonal polarizations, 6 represents a random ro- 
tation of the polarizations, whereas tpo and ipi denote the relative phases between the LO 
and the two modes, respectively. 

III. NUMERICAL RESULTS FOR TWO MODES 

In this section we present some Monte-Carlo simulations in order to estimate the working 
experimental conditions for performing the single-LO tomography on two-mode fields. We 
focus our attention on the twin-beam state, usually generated by spontaneous parametric 
downconversion, namely 

oo 

\^) = il-m'^J2^'ln)M,. (24) 

n=0 

For the simulations we use the following homodyne probability distribution that is derived 
in the Appendix: 

^27rA2(^,V'o,^i) V 2A2(^,V;o,V^i); 

where the variance A'i{6, tpo, tpi) is given by 

A2(n., ,. . l + |eP + 2|e|sin2gcos(^o + ^i-argO l-y . . 
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In the case of two radiation modes having the same frequency but orthogonal polarizations, 
Eq. (|25| ) gives the theoretical probability of outcome x for the homodyne measurement at 
a polarization angle 9, ipQ and ipi denoting the relative phases between the LO and the two 
modes, respectively. 

We study the tomographic measurement of the joint photon-number probability dis- 
tribution and the probability distribution for the total number of photons with use of the 
estimators in Eqs. (^) and (0), respectively. Moreover, we reconstruct the matrix elements 

Cn,m^a{m\b{m\^){^\n)a\n)b (27) 



that reveal the coherence of the twin-beam state by using the estimator in Eq. (|13D . For the 
twin-beam state in Eq. (p^, one should have 

Cn,m = il-\^\')rr- (28) 

The estimators have been numerically evaluated by applying the Gauss method for calcu- 
lating the integral in Eq. (|l^), which results in a fast and sufficiently precise algorithm with 
use of just 150 evaluation points. Notice that in the present case there is no convenience in 
using the factorization formula given in Ref. 0, as in that case an integral of a product of 
functions is needed. 

In Fig. |I] a Monte-Carlo simulation of the joint photon-number probability distribution 
is presented. The simulated values compare very well with the theoretical ones. We have 
done a careful analysis of the statistical errors for various twin-beam states by constructing 
histograms of deviations of the results from different simulated experiments from the the- 
oretical ones. In comparison to the customary two-LO tomography |T^, where for r/ = 1 
the statistical errors saturate for increasingly large n and m, here we have statistical errors 
that are slowly increasing versus n and m. This is due to the fact that the range of the 
estimator in Eq. (^) increases versus n and m. Overall we find that for any given quantum 
efficiency the statistical errors are generally slightly larger than those obtained with the 
two-LO method. The convenience of using a single LO then comes with its own price tag. 

By using the estimator in Eq. (|19D we have also constructed the probability distribution 
for the total number of photons N of the twin-beam state with unity (Fig. ^ as well as 
less-than-unity (Fig. ^ quantum efficiencies. Notice the dramatic increase of error bars 
versus N and for smaller rj. Finally, in Fig. ^ we report the results of the tomographic 
measurement of the matrix elements Cn,m defined in Eq. (|27|) . Because the reconstructed 
Cn,m is close to the theoretically expected value in Eq. (^), these reveal the purity of the 
twin-beam state, which cannot be inferred from the thermal diagonal distribution of Fig. |l|. 

IV. AN APPLICATION TO THE GHZ STATE 

Multimode homodyne tomography allows one to verify the generation of multimode 
states that are of interest in studies of the foundations of quantum mechanics. An example 
is the Greenberger-Horne-Zeilinger (GHZ) state |T^, which is a 6-mode state given by 



|GHZ) = ^ {\lao Ibo ICo) - \lae IK ICe)) , (29) 
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wherein o, e denote a couple of orthogonal polarizations; a, b, c pertain to electromagnetic 
modes with different wavevectors and/or frequencies; and the notation 1 1,1,1) represents 
the tensor product of three single-photon Fock states. The GHZ state is very interesting as 
it leads to correlations between three particles that are in contradiction with the Einstein- 
Podolsky- Rosen idea of "elements of reality" |18[. We note here that no experiment has yet 



succeeded in realizing the GHZ state. 

A tomographic measurement of the state in Eq. ( ppj ) can be suitably performed by 
varying randomly the phases and polarizations of the pairs of modes a^e, ^o.e, and Co^e, and 
then collecting homodyne outcomes by using three different LO's. The need of using three 



separate LO's in the present case is that in the actual experimental arrangement |T^ the 
three beams corresponding to modes Ooe, &o,e, and Co,e come with different wave- vectors 
and thus are spatially separated. Hence, such an experimental arrangement here gives the 
opportunity of using a combination of the present multimode method and the usual many- 
LO method based on the product of single-mode estimators. 

A simple tomographic check of the GHZ-state production consists of measuring the 
expectation value on the projector \(p){(p\, where 

\(P) = ^{\laolKlCo) + e'^\laelKlCe)) , (30) 

and comparing the result with the theoretical value, namely, 

C(0) = |(0|GHZ)p = i(l-cos0) . (31) 

Notice that for (f) = it the function C{(f)) represents the fidelity of the GHZ-state production. 
In addition, the same set of homodyne data allows one to recover the whole interference 
profile in Eq. (|3TD for all values of 0. 

In Fig. ^ we report the results of a Monte-Carlo simulation of the tomographic measure- 
ment of C(0) in Eq. (|3TD. We used 5 x 10^ homodyne data samples and assumed a quantum 
efficiency r] = 85%. For these parameters, the simulated C(0) compares very well with the 
theoretical one. 



V. CONCLUSIONS 

We have presented a generalization of the quantum homodyne tomography method to 
many modes of the radiation field that requires the use of only a single LO. By varying 
suitable random parameters the LO scans over all the linear combinations of the field modes. 
We have also provided a general method to obtain the "unbiased estimator" for a generic 
multimode operator. The quantum expectation value of such an operator can be evaluated 
for any unknown state of the radiation field by averaging the estimator over the homodyne 
outcomes that are collected by using a single LO. The estimators for some observables, 
such as the matrix elements of the multimode density operator and the total number of 
photons, have been explicitly evaluated. For the two-mode case we derived the estimator for 
the four- dimensional Q-function and the moments generating function of the total number 
of photons. By means of Monte-Carlo simulations we have analyzed in detail the case 



of the twin-beam state, namely the two-mode state produced by nondegenerate parametric 
amphfication (spontaneous downconversion) . For quantum efficiency of homodyne detection 
in the 80-90% range and with number of data samples of order 10^-10^, the simulations 
show that measurements of the joint photon-number probability, the distribution of the 
total number of photons, and the density-matrix elements are experimentally feasible. 

We have also shown an application of the method of multimode homodyne tomography 
to the measurement of the radiation field prepared in the Greenberger-Horne-Zeilinger state. 
The results of our simulations suggest that with a number of homodyne data samples around 
10^ and a homodyne detection efficiency of 85% the method would allow the reconstruction 
of such an interesting state of the radiation field with relatively small statistical errors. 

APPENDIX 

In this appendix we derive the theoretical probability distribution p,y(x; 9,ipQ,ipi) of the 
twin-beam state 

oo 

\^) = ^(x)|0)a|0), = (1 - m"^Y.^^n)M), , (32) 

n=0 

where S{x) = exp(xa"''6"'' — x*a6) and ^ = e*'^''^-^ tanh|x|. For unity quantum efficiency, the 
probability density p(x; 6, t/'q, V'l) is defined as follows: 

p(x;0,7/;o,^i) = Tr[f/t|a;),a(a;|®l,f/|^)(^|] 

= a(0|5(0| S\x) U^ \x)aa{x\ ® l,U S{x) |0)jO), , (33) 

where \x)a is the eigenvector of the quadrature x = |(a"'" + a) with eigenvalue x and U is the 
unitary operator achieving the mode transformation 

t (o\ jj _ ( e-''^° cos e e-'^' sin 9 \ fa\ 
VV \-e'^'sm9 e^'^^cos^ J \b) ' ^ ^ 

In the case of two radiation modes having the same frequency but orthogonal polarizations — 
the case of Type II phase-matched parametric amplifier — Eq. (^) gives the theoretical 
probability of outcome x for the homodyne measurement at a polarization angle 6 with 
respect to the polarization of the a mode, and with ipo and ipi denoting the relative phases 
between the LO and the two modes, respectively. By using the Dirac-5 representation of 
the X-quadrature projector 

— exp[iA(X — x)] , (35) 

Eq. (P^D can be rewritten as follows: 



^a(0|,(0|5t(x)f/te^^(^-^')f/5(x)|0)jO),= / ^e-^^x (36) 
a(0|b(0| exp U- [(e-^'^Vcos^ + e^'^V* sin^)a + {e'^^v* cosO + e-'"^' iismO)}) + H.c] | \Q)a\Q)b , 
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where we have used Eq. (0) and the transformation 



with /i = cosh|x| and v = e*'^''^-^ sinh|x|. Upon defining 

KC = e-'^'>fi cos e + e'^' u* sin 6 , 

KD = e'^° V* cos 6 + e"^'^^ /x sin ^ , (38) 

where iT G M and C, D e C, with |C|2 + \D\'^ = I one has 

ir^ =/i2_^ |z/p + 2/i|i/|sin26'cos(V'o + V'i -argz/) . (39) 

Now, since the unitary transformation 

(40) 

has no effect on the vacuum state, Eq. ([36| ) leads to the following Gaussian distribution: 

^ e-*"^(OUO| exp jzir- [(C a + Z^ 6) + H.c] | |0)jO)t 
— e-^^\(0| exp Uk^ [a + at] j |0), = - U(0|x/JOa|' 
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v/27rA2(^,V;o,V^i)''''^V 2A2(^, ^^0,^1) J ' ^^^^ 

where the variance A'^{9,il)o,%pi) is given by 

.2(n. -^ ^' l + |eP + 2|e|sin2gcos(V;o + ^i-argO ,,.. 

A (^, V^o, V^i) = — = 4(i_|^|2) • (42) 

Taking into account the Gaussian convolution that results from less-than-unity quantum 
efficiency, the variance just increases as A2(6', ipQ, ipi) -^ A^(^, ipQ, ipi) = ^'^{6, ipQ, ipi) + ^j^. 
Notice that the probability distribution in Eq. (^T|) corresponds to a squeezed vacuum for 
6^ = I and iIjq + ipi — arg^ = or vr. 
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FIG. 1. Two-mode photon-number probability p{n,m) of the twin-beam state of parametric 
fluorescence in Eq. ( |2^ ) for average number of photons per beam n = |^p/(l — |,^p) = 5 obtained 
by a Monte-Carlo simulation of the probability in Eqs. (^) and ( P6|) with random parameters 
cos 20, tpo, and tjji, using the estimator in Eq. (p^. On the left we have quantum efficiency r/ = 1 
and 10^ data samples were used in the reconstruction. On the right quantum efficiency ry = 0.9, 
and 5 x 10^ data samples were used. The theoretical values of off-diagonal p{n, m) are zero; for a 
comparison between theoretical and experimental diagonal p(n, n) probabilities and their relative 
statistical errors, see analogous experiments in Figs. |2| and ^. 
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FIG. 2. Probability distribution for the total number of photons of the twin-beam state 
in Eq. (|2^ ) for average number of photons per beam n = 2 obtained using the estimator in 
Eq. (p!9|). The oscillation of the total photon- number probability due to the perfect correlation of 
the twin-beam state has been reconstructed by simulating 10^ data samples with quantum effi- 
ciency rj = 1. The theoretical probability (thick solid line) is superimposed onto the result of the 
Monte-Carlo experiment; the latter is shown by the thin solid line with the statistical errors in 
gray shade. 
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FIG. 3. Similar to Fig. |^, but for quantum efficiency rj = 0.9 and 10^ data samples (on the 
left), and r] = 0.8 and 2 x lO'^ data samples (on the right). Notice the dramatic increase of error 
bars (in gray shade) versus N and for smaller rj. 
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FIG. 4. Tomographic reconstruction of the matrix elements Cn,m = aii^^lbi^^) {^Iit-) aln) b of 
the twin-beam state of parametric fluorescence in Eq. (|2^ ) for average number of photons per beam 
n = 2, obtained using the estimator in Eq. (|l3|). On the left we used 10^ simulated data samples 
and quantum efficiency rj = 0.9; on the right 3 x 10® data samples and rj = 0.8. The coherence of 
the twin-beam state is easily recognized as Cn,m varies little for n + m = constant [^ in Eq. ( p^ ) has 
been chosen real]. For a typical comparison between theoretical and experimental matrix elements 
and their relative statistical errors, see experiments in Figs. |2| and |^. 
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FIG. 5. Tomographic measurement of the overlap C((/>) between the GHZ state in Eq. (|29| ) and 
the state \(f)) in Eq. pO| ) with varying phase (/>. The value for cp = tt represents the fidelity between 
the experimental state and the theoretical one. Here a Monte-Carlo simulation with A^ = 2.5 x 10^ 
data samples and quantum efficiency t] = 0.85. The bars represent the statistical error, whereas 
the solid line is the theoretical value of C{(j)). All points are obtained from the same data samples 
(which causes the evident correlation between the statistical deviations). 
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